#ifndef GTL_MATH_HEADER_GUARD
#define GTL_MATH_HEADER_GUARD

//#include <cmath>
#include <math.h>
#include "gtlDebug.h"
#include "gtlCompilerSpecific.h"
#ifndef __APPLE__
#include "external/pstdint.h"
#endif

/// TODO: Endianness

namespace GTL
{
  namespace Math
  {
    ///************************************************************************/
    ///* Epsilons                                                             */
    ///************************************************************************/

    /// Function: EpsilonValue<T>()
    /// If type `T` is real, returns the smallest value larger than 0 for the type.
    /// If type `T` is integer, returns 0.
    template <typename T> T EpsilonValue();
    template <> inline float EpsilonValue<float>() { return 1.192092896e-07F; }
    template <> inline double EpsilonValue<double>() { return 2.2204460492503131e-016; }
    template <> inline int EpsilonValue<int>() { return 0; }

		inline bool CompareFloat(float a1, float a2, int maxUlps = 4)
		{
			union f2i { int i; float f; };

			f2i aInt; aInt.f = a1;
			if (aInt.i < 0)
				aInt.i = 0x80000000 - aInt.i;
			f2i bInt; bInt.f = a2;
			if (bInt.i < 0)
				bInt.i = 0x80000000 - bInt.i;
			return abs(aInt.i - bInt.i) <= maxUlps;
		}


    ///************************************************************************/
    ///* Square roots                                                         */
    ///************************************************************************/

    inline float SquareRoot(float v) { return sqrtf(v); }
    inline double SquareRoot(double v) { return sqrt(v); }
    inline int SquareRoot(int v) { return (int)sqrtf((float)v); }
    inline float SquareRootSSE(float x);
    template <typename T> T SquareRoot(const T& v) { return T::SquareRoot(v); }

    inline float InverseSquareRoot(float v) { return 1.0f/sqrtf(v); }
    inline double InverseSquareRoot(double v) { return 1.0/sqrt(v); }
    inline int InverseSquareRoot(int v) { return (int)(1.0f/sqrtf((float)v)); }
    inline float InverseSquareRootSSE(float x);

    ///************************************************************************/
    ///* Number manipulations                                                 */
    ///************************************************************************/

    template <typename T> T AbsoluteValue(T x) { return x>((T)0)?x:-x; }
    template <typename T> T Minimum(T a, T b) { return (a)<(b)?(a):(b); }
    template <typename T> T Maximum(T a, T b) { return (a)>(b)?(a):(b); }

    template <typename T> T ClampMinMaxValue(const T& minval, const T& maxval, const T& val) { return Maximum<T>(minval, Minimum<T>(maxval, val)); }

    /// Function: NextPowerOf2
    /// Returns the smallest power of 2 larger or equal to `v`
    inline uint32_t NextPowerOf2(uint32_t v) {
      v--; v |= v >> 1; v |= v >> 2; v |= v >> 4; v |= v >> 8; v |= v >> 16; v++; return v;
    }

		template <typename T>
		T SmoothStepMinMaxValue(const T& minval, const T& maxval, const T& val) {
			T clamped = ClampMinMaxValue((T)0, (T)1, (val - minval)/(maxval - minval));
			return clamped*clamped*(((T)3)-((T)2)*clamped);
		}

		template <typename T>
		T SmoothInterpolateStartEndAmount(const T& start, const T& end, const T& amount)
		{
			T clamped = ClampMinMaxValue((T)0, (T)1, amount);
			return start + ((end - start) * (clamped*clamped*(((T)3)-((T)2)*clamped)));
		}

    template <typename AB_TYPE, typename T_TYPE>
    AB_TYPE LinearInterpolateABT(const AB_TYPE& a, const AB_TYPE& b, const T_TYPE& t) {
      return a + (b - a) * ClampMinMaxValue((T_TYPE)0,(T_TYPE)1,t);
    }

		template <typename VALUE_TYPE, typename T_TYPE>
		VALUE_TYPE BilinearInterpolate(const VALUE_TYPE* fourCornersClockwise, const T_TYPE& x, const T_TYPE& y);

		template <typename T>
		inline T LinearToTriangle(const T& value, const T max_value, const T out_min, const T out_max)
		{
			return out_min + (1.0f - AbsoluteValue(((value/max_value) - 0.5f)*2.0f)) * (out_max-out_min);
		}

    template <typename T> T Floor(const T& v);
    template <> inline float Floor(const float& v) { return floorf(v); }
    template <> inline double Floor(const double& v) { return floor(v); }

    template <typename T> T Ceiling(const T& v);
    template <> inline float Ceiling(const float& v) { return ceilf(v); }
    template <> inline double Ceiling(const double& v) { return ceil(v); }

    template <typename T> T Sine(const T& v);
    template <> inline float Sine(const float& v) { return sinf(v); }
    template <> inline double Sine(const double& v) { return sin(v); }

    template <typename T> T Cosine(const T& v);
    template <> inline float Cosine(const float& v) { return cosf(v); }
    template <> inline double Cosine(const double& v) { return cos(v); }

		template <typename T> T ArcTan(const T& x, const T& y);
		template <> inline float ArcTan(const float& x, const float& y) { return atan2f(y, x); }
		template <> inline double ArcTan(const double& x, const double& y) { return atan2(y, x); }

    template <typename T> bool IsFinite(const T& v);
    template <> inline bool IsFinite(const float& v) { return v==v; } /// this ONLY works if we're using IEEE floats
    template <> inline bool IsFinite(const double& v) { return v==v; } /// this ONLY works if we're using IEEE floats
    template <> inline bool IsFinite(const int&) { return true; }

    template <typename T>
    T BendLinearCurveTowards(const T& value, const T& towards)
    {
      T exponent = log( towards ) * (T)-1.4427; /// (-1.4427 = 1 / log(0.5))
      return pow( value, exponent );
    }

    template <typename T> T BendLinearCurveTransitionTowards(const T& value, const T& towards);
    template <typename T> T SmoothBellCurve(const T& value);

    inline void SineCosine( float radians, float * GTL_RESTRICT_POINTER sine, float * GTL_RESTRICT_POINTER cosine );

    /// Does basically this to a number: 
    ///   0              m
		/// 0 #...............
		///   ........#.......
		///   ....#...........
		///   ............#...
		///   ..#.............
		///   ..........#.....
		///   ......#.........
		///   ..............#.
		///   .#..............
		///   .........#......
		///   .....#..........
		///   .............#..
		///   ...#............
		///   ...........#....
		///   .......#........
		/// v ...............#
    /// where m is maximum and v is value.
		/// `maximum` should be a power of two
    inline int32_t InsideOut(int32_t maximum, int32_t value);

    inline int _GreatestCommonDivisor (int i1, int i2)
    {
      if (i1 > i2)
      {
        if (i2 == 0)
          return (i1);
        return _GreatestCommonDivisor (i2, i1 % i2);
      }
      else
      {
        if (i1 == 0)
          return (i2);
        return _GreatestCommonDivisor (i1, i2 % i1);
      }
    }

		inline int GreatestCommonDivisor_UNTESTED(int i1, int i2)
		{
			for (;;)
			{
				if (i1 > i2)
				{
					if (i2 == 0) return (i1);
					int temp = i1%i2;
					i1 = i2;
					i2 = temp;
				}
				else
				{
					if (i1 == 0) return (i2);
					i2 = i2%i1;
				}
			}
		}

		inline bool IsPrime(unsigned int m)
		{
			unsigned int i, j;
			if (m < 2) return false;
			if (m == 2) return true;
			if (!(m & 1)) return false;
			if (m % 3 == 0) return (m == 3);
			for (i=5; (j=i*i), j <= m && j > i; i += 6) {
				if (m %   i   == 0) return false;
				if (m % (i+2) == 0) return false;
			}
			return true;
		}

    template <typename T>
    struct Constants { };

    template <>
    struct Constants<float>
    {
      static float Pi() { return 3.14159265358979f; }
      static float PiHalved() { return 3.14159265358979f/2.0f; }
      static float PiDoubled() { return 3.14159265358979f * 2.0f; }
      static float DegreesToRadians() { return 3.14159265358979f/180.0f; }
      static float RadiansToDegrees() { return 180.0f/3.14159265358979f; }

      static float FullCircleInRadians() { return PiDoubled(); }
      static float HalfCircleInRadians() { return Pi(); }
    };

    template <typename T, size_t mGranularity>
    class SinCosTable
    {
    public:
      SinCosTable()
      {
        mSinCosTable = 0;
      }
      ~SinCosTable()
      {
        delete[] mSinCosTable;
      }
      void Create()
      {
        if (mSinCosTable == 0)
        {
          mSinCosTable = new T[mGranularity];
          mGranularityInInverseRadians = ((T)mGranularity)/Constants<T>::FullCircleInRadians();
          for (size_t table_element=0; table_element<mGranularity; ++table_element)
          {
            mSinCosTable[table_element] = Sine<T>((T)table_element / mGranularityInInverseRadians);
          }
        }
      }
      T Sin(const T& radians)
      {
        Create();

        return mSinCosTable[((size_t)(radians*mGranularityInInverseRadians))%mGranularity];
      }
      T SinInterpolate(const T& radians)
      {
        Create();

        size_t fr = (size_t)Floor(radians*mGranularityInInverseRadians);
        size_t cr = (size_t)Ceiling(radians*mGranularityInInverseRadians);
        return LinearInterpolateABT(
          mSinCosTable[fr%mGranularity], 
          mSinCosTable[cr%mGranularity],
          radians-(T)fr);
      }

      T Cos(const T& radians) { return Sin(Constants<T>::PiHalved() - radians); }
      T CosInterpolate(const T& radians) { return SinInterpolate(Constants<T>::PiHalved() - radians); }

      T Tan(const T& radians) { return Sin(radians)/Cos(radians); }
      T TanInterpolate(const T& radians) { return SinInterpolate(radians)/CosInterpolate(radians); }

    private:
      T* mSinCosTable;
      T mGranularityInInverseRadians;
      //size_t mGranularity;
    };

    typedef SinCosTable<float, 6283> SinCosTableFloat;

    ///************************************************************************/
    ///* Implementations                                                      */
    ///************************************************************************/

    template <typename T>
    T BendLinearCurveTransitionTowards(const T& value, const T& towards)
    {
      if( value < (T)0.5 )
        return (T)0.5f * BendLinearCurveTowards(2*value, 1-towards);
      else
        return (T)1 - (T)0.5f * BendLinearCurveTowards(2 - 2*value, 1-towards);
    }

    inline void SineCosine( float radians, float * GTL_RESTRICT_POINTER sine, float * GTL_RESTRICT_POINTER cosine )
    {
#if defined( GTL_COMPILER_VISUAL_STUDIO )
      _asm
      {
        fld		DWORD PTR [radians]
        fsincos

          mov edx, DWORD PTR [cosine]
        mov eax, DWORD PTR [sine]

        fstp DWORD PTR [edx]
        fstp DWORD PTR [eax]
      }
#elif defined( GTL_COMPILER_GCC )
      register double __cosr, __sinr;
      __asm __volatile__ ("fsincos" : "=t" (__cosr), "=u" (__sinr) : "0" (radians));

      *sine = __sinr;
      *cosine = __cosr;
#else
      *sine = sinf(radians);
      *cosine = cosf(radians);
#endif
    }

    inline int32_t InsideOut(int32_t maximum, int32_t value)
    {
      int32_t m = maximum;
      int32_t result = 0;

      for (int32_t k=1; k<maximum; k<<=1)
      {
        if (value<<1 >= m)
        {
          result += k;
          value -= (m+1)>>1;
          m >>= 1;
        } else m = (m+1)>>1;
      }

      return result;
    }

    inline float SquareRootSSE(float x)
    {
      float   root = 0.f;
#ifdef GTL_COMPILER_VISUAL_STUDIO
      _asm
      {
        sqrtss          xmm0, x
          movss           root, xmm0
      }
#elif defined(GTL_COMPILER_GCC)
      __asm__ __volatile__(
        "movss %1,%%xmm2\n"
        "sqrtss %%xmm2,%%xmm1\n"
        "movss %%xmm1,%0"
        : "=m" (root)
        : "m" (x)
        );
#endif
      return root;
    }

    template <typename T>
    T SmoothBellCurve(const T& value)
    {
      return (1 - Cosine( value * Constants<T>::Pi() )) * (T)0.5f;
    }

		template <typename VALUE_TYPE, typename T_TYPE>
		VALUE_TYPE BilinearInterpolate(const VALUE_TYPE* fourCornersClockwise, const T_TYPE& x, const T_TYPE& y)
		{
			T_TYPE omx = (T_TYPE)1 - x;
			T_TYPE omy = (T_TYPE)1 - y;
			return fourCornersClockwise[0] * omx * omy + 
				fourCornersClockwise[1] * x * omy +
				fourCornersClockwise[3] * omx * y +
				fourCornersClockwise[2] * x * y;
		}
  }
}

#endif